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LONG-TERM  GOALS 

The  long-term  technical  goal  is  to  design,  develop  and  test  the  next  generation,  primitive  equation 
ocean  model  for  high-resolution  scientific  (ROMS:  Regional  Ocean  Modeling  System)  and  operational 
(TOMS:  Terrain- following  Ocean  Modeling  System)  applications.  This  project  will  improve  the  ocean 
modeling  capabilities  of  the  U.S.  Navy  for  relocatable,  coastal,  coupled  atmosphere-ocean  forecasting 
applications.  It  will  also  benefit  the  ocean  modeling  community  at  large  by  providing  the  current 
state-of-the-art  knowledge  in  physics,  numerical  schemes,  and  computational  technology. 

OBJECTIVES 

The  main  objective  is  to  produce  a  tested  expert  ocean-modeling  framework  for  scientific  and 
operational  applications  over  a  wide  range  of  spatial  (coastal  to  basin)  and  temporal  (days  to  seasons) 
scales.  The  primary  focus  is  to  implement  the  most  robust  set  of  options  and  algorithms  for  relocatable 
coastal  forecasting  systems  nested  within  basin-scale  operational  models  for  the  Navy.  The  system 
includes  some  of  the  analysis  and  prediction  tools  that  are  available  in  Numerical  Weather  Prediction 
(NWP),  such  as:  4-dimensional  variational  data  assimilation  (4D-Var),  ensemble  prediction,  adaptive 
sampling,  and  circulation  stability  and  sensitivity  analysis. 

APPROACH 

The  structure  of  TOMS  is  based  on  ROMS  because  of  its  accurate  and  efficient  numerical  algorithms, 
tangent  linear  and  adjoint  models,  variational  data  assimilation,  modular  coding  and  explicit  parallel 
structure  conformal  to  modem  computer  architectures  (both  cache-coherent  shared-memory  and 
distributed  cluster  technologies).  Currently,  both  ROMS  and  TOMS  are  identical  and  continue 
improving  and  evolving.  ROMS  remains  as  the  scientific  community  model  while  TOMS  becomes  the 
operational  community  model. 

ROMS/TOMS  is  a  three-dimensional,  free-surface,  terrain-following  ocean  model  that  solves  the 
Reynolds-averaged  Navier-Stokes  equations  using  the  hydrostatic  vertical  momentum  balance  and 
Boussinesq  approximation  (Haidvogel  et  at.  2000,  2008;  Shchepetkin  and  McWilliams,  2005,  2009). 
The  governing  dynamical  equations  are  discretized  on  a  vertical  coordinate  that  depends  on  the  local 
water  depth.  The  horizontal  coordinates  are  orthogonal  and  curvilinear  allowing  Cartesian,  spherical, 
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and  polar  spatial  discretization  on  an  Arakawa  C-grid.  Its  dynamical  kernel  includes  accurate  and 
efficient  algorithms  for  time-stepping,  advection,  pressure  gradient  (Shchepetkin  and  McWilliams 
2003,  2005),  several  subgridscale  parameterizations  (Durski  et  al.,  2004;  Warner  et  al.,  2005)  to 
represent  small-scale  turbulent  processes  at  the  dissipation  level,  and  various  bottom  boundary  layer 
formulations  to  detennine  the  stress  exerted  on  the  flow  by  the  bottom. 

Several  adjoint-based  algorithms  exist  to  explore  the  factors  that  limit  the  predictability  of  the 
circulation  in  regional  applications  for  a  variety  of  dynamical  regimes  (Moore  et  al,  2004,  2009). 
These  algorithms  use  the  ideas  of  Generalized  Stability  Theory  (GST)  in  order  to  identify  the  most 
unstable  directions  of  state-space  in  which  errors  and  uncertainties  are  likely  to  grow.  The  resulting 
singular  vectors  can  be  used  to  construct  ensembles  of  forecasts  by  perturbing  initial  and  boundary 
conditions  (optimal  perturbations)  and/or  perturbing  surface  forcing  (stochastic  optimals).  Perturbing 
the  system  along  the  most  unstable  directions  to  the  state-space  yields  information  about  the  first 
(ensemble  mean)  and  second  (ensemble  spread)  moments  of  the  probability  density  function.  Given  an 
appropriate  forecast  skill  measure,  the  circulation  is  predictable  if  low  spread  and  unpredictable  if 
large  spread. 

ROMS/TOMS  uniquely  supports  three  different  4D-Var  data  assimilation  methodologies  (Moore  et 
al.,  2010a,  b):  a  primal  form  of  the  incremental  strong  constraint  4D-Var  (I4D-Var),  a  strong/weak 
constraint  dual  form  of  4D-Var  based  on  the  Physical-space  Statistical  Analysis  System  (4D-PSAS), 
and  a  strong/weak  constraint  dual  fonn  of  4D-Var  based  on  the  indirect  representer  method  (R4D- 
Var).  In  the  dual  formulations,  the  search  for  the  best  ocean  circulation  estimate  is  in  the  subspace 
spanned  only  by  the  observations,  as  opposed  to  the  full  space  spanned  by  the  model  as  in  the  primal 
formulation.  Although  the  primal  and  dual  formulations  yield  identical  estimates  of  the  ocean 
circulation  for  the  same  a  priori  assumptions,  there  are  practical  advantages  and  disadvantages  to  both 
approaches  (Moore  et  al.,  2010a,  b,  c).  To  our  knowledge,  ROMS/TOMS  is  the  only  open-source, 
ocean  community-modeling  framework  supporting  all  these  variational  data  assimilation  methods  and 
other  sophisticated  adjoint-based  algorithms. 

There  are  several  biogeochemical  models  available  in  ROMS.  In  order  of  increasing  ecological 
complexity  these  include  three  NPZD-type  models  (Franks  et  al,  1986;  Powell  et  al.,  2006;  Fiechter  et 
al,  2009),  a  nitrogen-based  ecosystem  model  (Fennel  et  al,  2006,  2008),  a  Nemuro-type  lower  level 
ecosystem  model  (Kishi  et  al,  2007),  and  a  bio-optical  model  (Bissett  et  al,  1999). 

ROMS  includes  a  sediment-transport  model  with  an  unlimited  number  of  user-defined  cohesive  (mud) 
and  non-cohesive  (sand)  sediment  classes  (Warner  et  al,  2008).  Each  class  has  attributes  of  grain 
diameter,  density,  settling  velocity,  critical  stress  threshold  for  erosion,  and  erodibility  constant.  A 
multi-level  bed  framework  tracks  the  distribution  of  every  size  class  in  each  layer  and  stores  bulk 
properties  including  layer  thickness,  porosity,  and  mass,  allowing  the  computation  of  bed  morphology 
and  stratigraphy.  Also  tracked  are  bed-surface  properties  like  active-layer  thickness,  ripple  geometry, 
and  bed  roughness.  Bedload  transport  is  calculated  for  mobile  sediment  classes  in  the  top  layer. 

ROMS  is  a  very  modem  and  modular  code  written  on  F90/F95.  It  uses  C-preprocessing  to  activate  the 
various  physical  and  numerical  options.  The  parallel  framework  is  coarse-grained  with  both  shared- 
memory  (OpenMP)  and  distributed-memory  (MPI)  paradigms  coexisting  in  the  same  code.  Because  of 
its  construction,  the  parallelization  of  the  adjoint  is  only  available  for  MPI.  Several  coding  standards 
have  been  established  to  facilitate  model  readability,  maintenance,  and  portability.  All  the  state  model 
variables  are  dynamically  allocated  and  passed  as  arguments  to  the  computational  routines  via 
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dereferenced  pointer  structures.  All  private  arrays  are  automatic;  their  size  is  determined  when  the 
procedure  is  entered.  This  code  structure  facilitates  computations  over  nested  grids.  There  are  three 
types  of  nesting  capabilities  in  ROMS:  (i)  refinement  grids  which  provide  increased  resolution  (3:1  or 
5:1)  in  a  specific  region;  (ii)  mosaics  which  connect  several  grids  along  their  edges,  and  (iii)  composite 
grids  which  allow  overlap  regions  of  aligned  and  non-aligned  grids. 

WORK  COMPLETED 

After  a  complete  overhaul  of  the  ROMS/TOMS  4D-Var  data  assimilation  capabilities  in  FY09,  we 
concentrated  in  FY10  to  develop  and  test  additional  algorithms  for  4D-Var  observations  impact, 
observations  sensitivity,  array  modes,  and  posterior  error  covariance  analysis.  Some  of  these 
developments  required  the  adjoint  of  the  4D-Var  data  assimilation  system,  denoted  as  (4DVar)T. 
Currently,  (4DVar)T  is  only  available  for  the  dual  formulation  algorithms  4D-PSAS  and  R4D-Var 
because  the  adjoint  of  the  Lanczos-based,  conjugate  gradient  algorithm  used  in  the  minimization  is 
much  simpler.  This  capability  is  challenging  in  I4D-Var  due  to  the  complexity  of  the  I/O  intensive, 
Lanczos-based,  conjugate  gradient  and  preconditioning  algorithms  that  are  used  when  the 
minimization  is  carried  out  in  model  space.  We  are  planning  to  build  the  (I4D-Var)T  in  the  future. 

The  observation  impact  algorithm  can  be  used  to  quantify  the  contribution  of  each  observation  during  a 
4D-Var  analysis  to  a  specified  aspect  (scalar  function,  say  I)  of  the  ocean  circulation  (Langland  and 
Baker,  2004;  Gelaro  and  Zhu,  2009;  Tremolet,  2008;  Moore  et  al.,  2010c).  That  is,  it  identifies  the  part 
of  the  model  space  that  controls  /  and  that  is  activated  by  the  observations.  It  yields  the  actual 
contribution  of  each  observation  to  the  circulation  increment. 

The  observation  sensitivity  is  based  on  (4D-Var)T  and  quantifies  the  change  that  would  result  in  the 
circulation  estimate  as  result  of  changes  in  the  observations  or  observation  array  (Tremolet,  2008; 
Moore  et  al,  2010c).  It  is  a  very  useful  tool  for  efficient  generation  of  observation  system  experiments 
(OSEs),  observation  array  design,  and  adaptive  sampling.  It  also  can  be  used  to  predict  the  changes 
that  will  occur  in  I  in  the  event  of  a  platform  failure  or  a  change  in  the  observation  array  (Moore  et  al. , 
2010c). 

The  array  modes  of  the  stabilized  representer  matrix  can  be  used  to  determine  the  most  stable 
component  of  the  circulation  with  respect  to  changes  in  the  innovation  vector  (Bennett,  1985;  Moore  et 
al.,  2010c).  The  array  modes  are  independent  of  the  observation  values  and  depend  only  on  the 
observation  locations,  the  prior  covariances,  and  prior  circulation. 

Substantial  work  has  been  carried  out  in  the  nesting  capabilities  of  ROMS/TOMS.  Currently,  we  are 
migrating  this  development  from  our  research  branches  to  the  distributed,  official  version  of  the  model. 
The  kernel  of  the  nonlinear,  tangent  linear,  representer,  and  adjoint  models  has  been  updated  to  allow 
different  I-  and  J-  ranges  in  DO-loops  to  permit  operations  on  various  nested  grid  classes  (refinement, 
mosaics,  and  composite)  and  nesting  layers  (refinement  and  composite  grid  combinations).  This 
development  is  very  delicate  and  requires  a  lot  of  testing.  We  have  been  developing  this  capability  for 
couple  of  years  and  will  be  available  soon  to  the  user  community  and  released  as  ROMS  4.0. 

We  held  a  very  successful  workshop  at  the  Hawaii  Imin  International  Conference  Center,  University  of 
Hawaii  at  Manoa,  Honolulu,  Hawaii,  April  5-8,  2010  (http://www.myroms.org/2010workshop).  As  in 
the  past,  several  tutorials  were  offered  about  basic  and  advanced  ROMS  algorithms. 


3 


We  also  organized  a  4D-Var  data  assimilation  workshop  for  advanced  ROMS  users  at  the  Simularium, 
Baskin  School  of  Engineering,  University  of  California  Santa  Cruz,  Santa  Cruz,  CA,  July  12-16,  2010 
(http://www.myroms.org/4DVar_workshop_2010).  It  included  tutorial  style  lectures  on  4D-Var  data 
assimilation  methodology,  algorithms,  and  hands-on  exercises.  By  the  end  of  the  week,  users  were 
running  their  own  application.  This  was  our  first  workshop  of  this  kind  and  serves  as  a  prototype  for 
similar  workshops  in  the  future. 

RESULTS 

A  comprehensive  documentation  material  was  developed  for  the  4D-Var  data  assimilation  workshop 
held  at  University  of  California  Santa  Cruz  in  July  2010.  The  lectures,  tutorial,  and  exercise  are 
available  to  everyone  in  WikiROMS,  as  shown  in  Figure  1. 


•  Q  https://www.myroms.org/wiki/index.php/4DVar_Tutorial_Introductjon 
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4-Dimensional  Variational  (4D-Var)  Data  Assimilation  Tutorial 
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This  tutorial  shows  several  examples  of  ROMS  4D-Var  data  assimilation  algorithms  in  the  California  Current  System  (CCS).  1/3  degree  resolution  (  ?). 
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:  model  domain  is  shown  in  Fig.  1  and  has  open  boundaries  along  the  northern,  western,  and  southern  edges  of  the  model  domain. 
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Fig.  1:  Model  Bathymetry  with  37*N  Transect 

In  the  tutorial,  you  wi  perform  a  4D-Var  data  assimilation  cycle  that  spans  the  period  3-6  January,  2004.  The  4D-Var  control  vector  6z  is  comprised  of  increments  to  the  initial  conditions,  6x(tg),  surface  forcing,  6f(f),  and  open 
boundary  conditions,  6b(f).  The  prior  initial  conditions.  X|,(t0),  are  taken  from  the  sequence  of  4D-Var  experiments  described  by  -re  e:  ■ .  in  which  data  were  assimilated  every  7  days  during  the  period  July  2002- 
December  2004.  The  prior  surface  forcing,  fb(f),  takes  the  form  of  surface  wind  stress,  heat  flux,  and  a  freshwater  flux  computed  using  the  ROMS  buk  flux  formulation,  and  using  near  surface  air  data  from  COAMPS  (Dc  vie  ef  si 
§*■  2C09).  Clamped  open  boundary  conditions  are  imposed  on  (w,v)  and  tracers,  and  the  prior  boundary  conditions,  b|,(f),  are  taken  from  the  global  ECCO  product  (‘.Vursch  ar  c  Heimbach  2007).  The  free-surface  height  and  vertically 
integrated  velocity  components  are  subject  to  the  usual  Chapman  and  Flather  radiation  conditions  at  the  open  boundaries.  The  prior  surface  forcing  and  open  boundary  conditions  are  provided  dafly  and  linearly  interpolated  in  time. 
Similarly,  the  increments  5f(f)  and  6b(t)  are  also  computed  daily  and  linearly  interpolated  in  time, 
ait  mooTM  19:tl.  14  July  Sit.  KC*MMSMt>n>M 


Figure  1:  Screenshot  of  WikiRoms  4D-Var  data  assimilation  tutorials  showing  downloadable  PDF 
links  for  lectures,  tutorials,  and  exercises.  It  includes  pages  for  1)  Introduction,  2)  Error  Covariance 
Normalization,  3)  I4D-Var,  4)  I4D-Var  Observation  Impact,  5)  4D-PSAS,  6)  4D-PSAS  Observation 
Impact,  7)  4D-PSAS  Observation  Sensitivity,  8)  R4D-Var,  9)  R4D-Var  Observation  Impact,  10) 
R4D-Var  Observation  Sensitivity,  and  11)  Array  Modes. 
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A  realistic  ROMS  application  for  the  U.S.  west  coast  and  the  California  Current  System  (CCS)  is 
provided  in  the  tutorials  to  run  all  the  4D-Var  algorithms  introduced  during  the  lectures.  This 
configuration  is  referred  to  as  WC13.  It  has  10  km  horizontal  resolution  and  30  levels  in  the  vertical. 
While  30  km  resolution  is  inadequate  for  capturing  much  of  the  energetic  circulation  associated  with 
the  CCS,  WC13  captures  the  broad  scale  features  of  the  circulation  quite  well,  and  serves  as  a  very 
useful  and  efficient  illustrative  example  for  the  various  4D-Var  data  assimilation  algorithms  in  the 
tutorials  for  the  period  of  Jan  3-6,  2004. 

The  4D-Var  control  vector  8z  is  comprised  of  increments  to  the  initial  conditions,  8x(t0),  surface 
forcing,  8f(t),  and  open  lateral  boundary  conditions,  8b(t).  The  prior  initial  conditions,  xh(t0),  are 
taken  from  the  sequence  of  4D-Var  experiments  described  by  Moore  et  al.  (2010b)  in  which  data  were 
assimilated  every  7  days  during  the  period  July  2002-  December  2004.  The  prior  surface  forcing, 
f b(t),  takes  the  form  of  surface  wind  stress,  heat  flux,  and  a  freshwater  flux  computed  using  the  ROMS 
bulk  flux  fonnulation,  and  using  near  surface  air  data  from  COAMPS  (Doyle  et  al,  2009).  Clamped 
open  boundary  conditions  are  imposed  on  (w,v)  and  tracers,  and  the  prior  lateral  open  boundary 
conditions,  b b(t),  are  taken  from  the  global  ECCO  product  (Wunsch  and  Heimbach,  2007).  The  firee- 
surface  height  and  vertically  integrated  velocity  components  are  subject  to  the  usual  Chapman  and 
Flather  radiation  conditions  at  the  open  boundaries.  The  prior  surface  forcing  and  lateral  open 
boundary  conditions  are  provided  daily  and  linearly  interpolated  in  time.  Similarly,  the 
increments  8f(t)  and  8b(t)  are  also  computed  daily  and  linearly  interpolated  in  time. 

The  observations  assimilated  into  the  model  are  satellite  SST,  satellite  SSH  in  the  fonn  of  a  gridded 
product  from  Aviso,  and  hydrographic  observations  of  temperature  and  salinity  collected  from  Argo 
floats  and  during  the  GLOBEC/LTOP  and  CalCOFI  cruises  off  the  coast  of  Oregon  and  southern 
California,  respectively.  The  observation  locations  are  illustrated  in  Figure  2. 

The  majority  of  the  observations  available  for  assimilation  in  the  ocean  are  from  satellite  platforms  in 
the  fonn  of  SSH  and  SST.  The  analysis  of  Moore  et  al.  (2010b)  found  that  less  than  10%  of  the 
available  satellite  observations  provide  independent  information.  There  is  a  lot  of  observation 
redundancy  which  increases  with  model  horizontal  resolution. 

Figure  3,  shows  the  observation  impact  and  observation  sensitivity  for  a  3  day  R4D-Var  data 
assimilation  cycle  on  the  time-averaged  transport  across  37N  over  the  upper  500m,  denoted  by  Bin, 
and  given  by: 

/37W(x)=^i>rx, 

i=l 

where  h  is  a  vector  with  non-zero  elements  corresponding  to  the  velocity  grid  points  that  contribute  to 
the  transport  normal  to  the  37N  section  shown  in  Figure  2,  N  is  the  number  of  time  steps  during  the 
assimilation  interval,  and  X,  is  the  model  state  vector  at  time  i  At .  Figure  3a  shows  the  total  number  of 
observations  and  number  observations  per  platform  used  in  the  observation  impact  and  sensitivity 
analysis  of  the  transport  across  37N.  It  is  clear  that  the  SST  is  the  dominant  observation  platform 
followed  by  SSH.  The  analysis  (Fig.  3b)  of  the  alongshore  transport  in  this  upwelling  area  is 
influenced  primarily  by  uncertainties  in  the  initial  conditions,  8x(t0),  and  surface  forcing,  8f(t). 

Notice  that  the  impact  of  the  CTD  and  Argo  observations  on  the  Imm  is  significant  despite  the  small 
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number  of  these  observations  (Fig.  3c).  The  observation  sensitivity  analysis  (Fig.  3d)  indicates  that 
/37v  will  be  affected  considerably  by  any  changes  in  the  SSH  coverage  (Moore  et  al.,  2010c). 
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Figure  2:  Observations  available  for  4D-Var  data  assimilation  during  the  period  of  Jan  3-6,  2004 : 
(a)  Aviso  Sea  Surface  Height  (SSH)  anomaly,  (b)  Sea  Surface  Temperature  (SST)  daily  composite 
from  NOAA  PFEG  CoastWatch,  (c)  in  situ  temperature,  and  (d)  in  situ  salinity  from  Argo  floats, 

and  GLOBE C/LTOP  and  CalCOFI  cruises. 
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Figure  3:  4D-Var  observation  impact  and  sensitivity  for  the  transport  across  37N:  (a)  total  number 
of  observations  and  number  of  observations  per  platform  used  in  the  observation  impact  and 
observation  sensitivity  analyses,  (b)  total  transport  increment  A 1 37 nz  (Sv)  and  increments  associated 
with  initial  conditions,  surface  forcing,  and  lateral  boundary  conditions,  (c)  same  as  (b)  except  that 
the  bars  indicate  the  contributions  from  each  observing  platform,  and  (d)  same  as  (c)  except  but  the 

observation  sensitivity  using  (4D-Var)T . 
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IMPACT/APPLICATIONS 


This  project  will  provide  the  ocean  modeling  community  with  a  freely  accessible,  well  documented, 
open-source,  terrain-following,  ocean  model  for  regional  nowcasting  and  forecasting  that  includes 
advanced  data  assimilation,  ensemble  prediction,  and  analysis  tools  for  adaptive  sampling  and 
circulation  dynamics,  stability,  and  sensitivity. 

TRANSITIONS 

The  full  transition  of  ROMS/TOMS  to  the  operational  community  is  likely  to  occur  in  the  future. 
However,  the  ROMS/TOMS  algorithms  are  now  available  to  the  developers  and  scientific  and 
operational  communities  through  the  website  http://www.myroms.org/. 

RELATED  PROJECTS 

The  work  reported  here  is  related  to  other  already  funded  ONR  projects  using  ROMS.  In  particular, 
the  PI  (H.  Arango)  closely  collaborates  with  A.  Moore  (adjoint-based  algorithms)  at  University  of 
California,  Santa  Cruz,  A.  Miller  and  B.  Cornuelle  (ROMS  adjoint  and  variational  data  assimilation)  at 
Scripps  Institute  of  Oceanography,  E.  Di  Lorenzo  (Southern  California  predictability)  at  Georgia 
Institute  of  Oceanography,  and  J.  Wilkin  (Mid- Atlantic  Bight  variational  data  assimilation)  at  Rutgers 
University. 

The  PI  is  also  supported  by  the  following  grants: 

“Characterization  and  Modeling  of  the  Philippine  Archipelago  Dynamics  Using  the  ROMS  4DVAR 
Data  Assimilation  System”,  grant  number  N00014-04-1-0417,  http://www.myroms.org/philex. 
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